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Abstract 

We performed numerical simulations of blood flow in arteries with a 
variable stiffness and cross-section at rest using a finite volume method 
coupled with a hydrostatic reconstruction of the variables at the interface 
of each mesh cell. The method was then validated on examples taken from 
the literature. Asymptotic solutions were computed to highlight the effect 
of the viscous and viscoelastic source terms. Finally, the blood flow was 
computed in an artery where the cross-section at rest and the stiffness 
were varying. In each test case, the hydrostatic reconstruction showed 
good results where other simpler schemes did not, generating spurious 
oscillations and nonphysical velocities. 


1 Introduction 

In this work we are interested in modeling and simulating blood flow in arteries 
with varying stiffness and cross-section. The blood flow in the main arteries of 
the systemic network is governed by the 3D Navier-Stokes equations which can 
be complicated and time-consuming to solve numerically. Fortunately, using 
well-known hypothesis valid in the case of blood flow in arteries (long wave 
approximation D/\ « 1, axial symmetry dg = 0), this system of equations can 
be simplified and then integrated over the cross-section of the artery in order to 
obtain a ID hyperbolic system of equations, similar to the Saint-Venant system 
for shallow water flows. Details on the derivation of the model are presented in 
section (2) and can also be found in [23, 38]. Finally, we are left with a set of 
mass and momentum conservation equations with non dimensionless variables 
and parameters: 
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with A{x,t) = TTR{x,t)^ the cross-section area {R is the radius of the artery), 
Q{x, t) = A{x, t)u{x, t) the discharge, u{t, x) the mean flow velocity, p the blood 
density, C/ the friction coefficient and Ao = k\fAQ with k{x) the stiffness of 
the artery and Ao{x) = 'nRo{xY the cross-section at rest. 

The vast majority of arteries in the systemic network are tapered, meaning 
that the cross-section at rest Aq (x) varies throughout the length of the artery. 
Similarly, in the presence of arterial pathologies such as aneurysm or stenoses, 
the stiffness k {x) of the arterial wall can vary locally. As for shallow water 
equations with topography, the presence of tapper or variable stiffness in an 
artery modifies the blood flow, and both behaviors are accounted for in (1) 
through the source term A ( — 2'/Adxkj3 ) l\pkp. To numerically solve 

(1), it is necessary, among other things, to discretize this source term. A naive 
treatment of the topography gradients will most likely generate numerical os¬ 
cillations, therefore the use of the so-called well-balanced schemes is required to 
properly balance the fluxes and the source terms. In the following, we will focus 
on a specific well-balance method, called the hydrostatic reconstruction. 

We will first present the derivation of the model and its properties, then the 
numerical method and in particular the derivation of the well-balanced scheme 
applied to the case of blood flow in arteries. We will then validate our method 
on examples taken from the literature and verify asymptotic behaviors of the 
numerical solution. Finally, we will compute the blood flow in an artery with 
varying cross-section and stiffness. 


2 Derivation of the ID blood flow equations 


The ID model for blood flow equations is derived from the conservative form of 
the Navier-Stokes equations for an incompressible fluid with constant viscosity 

p.: 


dtp + \7pu = 0 

(2) 

dtpu V • {puu -\- pi + t) = 0, 

(3) 


where u is the velocity vector, p the density, supposed constant, p the pressure 
and r the stress tensor to be defined. Using the control volume of the Figure 
1, we integrate the Navier-Stokes equations over a volume V of cross-section A 
surrounded by a surface S (V = S U A) and of length dz. We define then the 
average velocity U and pressure P as 

{U,P} = 4 / {u,p}dA. 

^ JdA 

From the mass conservation equation (2) we have: 


{y pu)dV = / pu ■ ndS + 


pu 


ndA. 


IdV 


Ids 


IdA 


We then transform the volume integral using the Green (Divergence) theorem 
and writing the surface integral as S' U A. The surface element is dS = RdOdz 
and the two terms are written as 



• ndS = 27r 


Ur\RRdx 



^Rdx = 
at 




2 


V 



Figure 1: Control volume for integration (see text). 


and 


I pu ■ ndA = {pAU)i — {pAU )2 = j d{pAU) = f 

JdA J Jdx 

We retrieve therefore the first equation of our system 

dtA + dMU)=Q. 


dpAU 

dx 


For the conservation of momentum equation (3), the temporal term dtpu be¬ 
comes 


j dt{pu)dV = P [ dtudAdx = p j dt{UA)dx 

JdV ^dV J dx 

and the divergence term 

Jqv ^ ■ {puu -\- pi + t) = / {puu -\- pi t) ■ ndS + 

Jas 

/ {puu + pi + t) ■ ndA. 

JdA 

In the last two integrals the integration over the surface S is 
/ {puu + pi t) ■ ndS = / {pux + Trx)dS, 

Jas Jas 

where the term uudS tends to zero. Finally, the integration over the area A 
gives 


fg^{puu + pi + t) ■ ndA 


— [A{pU‘^ + P + Txx)]'i 

r dA{u^ + p/p) 

Idx 


= P 


dx. 


In terms of the cross-section A and the flow rate Q, we obtained the following 
system of equations: 


dtA + dxQ — 0 

dtQ + dx^ = -jdxP- h- 
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The viscous effects are contained in fy which is computed by the integration 
of the shear stress at the wall Tyx over the internal surface dS. Therefore, it 
depends on the exact flow condition. To close the mathematical problem we 
need a relation between the pressure P and the cross-section A, P = P{A), 
called the wall or state law. For fy = CfQjA and the state law P = Pq + 
k (x) l^/^{^/A{x,t) - {x)), which corresponds to the elastic response of 

the artery, we obtain the proposed system of equations. (1). 


3 Conservative hyperbolic system and steady states 

Considering an artery with a constant stiffness k and a variable cross-section 
at rest Aq{x), (1) reduces to the following system, similar to the shallow water 
equations with topography: 

i dtA dxQ = 0 

O' 

A 

As a reminder, the shallow water system is: 

( dth + dxq = 0 

I dtq + dx[^ + = gh (5o - Sf) , 

with h{x,t) the water height, q[x,t) = h(x,t)u{xA) the unit discharge, u{x,t) 
the mean flow velocity, g the constant of gravity, Sq = —dxZ the opposite of 
the slope, z the topography and Sf the friction term (which takes the form of 
Manning’s, Stickler’s, Chezy’s, ... empirical friction law). 




A3/2 = 


A 

pV^ 


dxAo - Cf 


(5) 


dtQ + dx 


3.1 Hyperbolic system 

The system (5) can be written using the following vectorial form: 

dtU + dxF{U) = S{U), (7) 


where U is the vector of the conservative variables, F{U) is the flux: 


U = 



A 


Zp^A?!'^ j 


( 8 ) 


and S{U) is the source term, taking into account the shape of the vessel at rest 
Aq (x) and the friction term 


S{U) 


pV^ 


dxAo-Cf^ 


(9) 


The analogous term for the shallow water equations is the topography source 
term. The gradient of the flux (8) can be written as the product of the Jacobian 
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matrix J{U) with the partial derivative of the vector of conservative variables 
U: 

d,F{U) = I kVA Q2 2Q ] .dj ^) = J{U).d,U. (10) 

V 2p^ A j ^ ^ ^ 

When the cross-section A > 0, the Jacobian matrix admits two different real 
eigenvalues, Ai and A 2 : 


A 

Ai--- 


I kVA ^ ^ Q 

-—= = u — c and A 2 = -r 
2pv^ A 


I k^/A 

2pv^ 


= u + c, 


( 11 ) 


with c the Moens-Korteweg wave propagation velocity (for the shallow water 
equations (6), c = \/gh). In this case, the system is said to be strictly hyperbolic, 
which is a generalization of the advection phenomenon [18, 45, 30]: a part of 
the information concerning the flow propagates at the velocity Ai and the other 
part at the velocity A 2 . For blood flow under physiological conditions, we have 
Ai > 0 and A 2 < 0, hence the flow is subcritical. 


3.2 Steady states 


Since the works of [3, 2] on the shallow water equations, it is well known that if a 
numerical scheme does not preserve steady states at the discrete level, spurious 
oscillations and artificial non zero velocities will be generated. The steady states 
for the system (5) are obtained when considering a stationary flow {i.e. there 
is no evolution in time) and are governed by the following equations: 


dxQ = 0 


( 12 ) 


with b = fc/(pv^) constant since we are considering an artery with a constant 
stiffness k. Neglecting the viscous friction effects (inviscid flow) by setting C/ = 
0, we obtain the conservation of the discharge and Bernoulli’s law for blood 
flow: 


Q = Qo 

^ -f bVA - by/Ao (x) = cst. 


(13) 


In the literature [8, 37, 44, 6], we can find well-balanced numerical methods able 
to preserve the following steady state: 


{ 9 = 90 

^ + h + z{x)=cst, 

which is the analogous of (13) in the case of the shallow water equations. How¬ 
ever, these methods are complicated to handle due to the occurrence of critical 
points when solving (13) or (14). Therefore we chose to focus on simpler steady 
states that we call the rest steady states or the ’’man at eternal rest” equilib¬ 
rium [13] by analogy with the ’’lake at rest” (introduced in [1]) or the hydrostatic 
equilibrium for the shallow water equations: 


f q = u = 0 

\ dxih + z (x)) = dxV = 0, 


(15) 
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where rj is the water level. In this case we have a hydrostatic balance between 
the hydrostatic pressure and the gravitational acceleration. By analogy, we have 
the following equilibrium for the blood flow in arteries: 

f Q = u = 0 

I dx (bVA - b^yAo (a;)^ = 0 . 

Numerical methods able to preserve at least the steady states (16) are said to be 
’’well-balanced” since the work of [19]. A wide panel of well-balanced methods 
has been developed for shallow water equations. Among others we can mention 
[29, 24, 39, 27, 16, 25, 1, 11, 36, 17, 4, 22, 5, 20]. In [13], we adapted the hydro¬ 
static reconstruction introduced in [1] to the system with constant stiffness (5). 


We will now present the hydrostatic reconstruction introduced in [1] adapted 
to the original system of equations (1) with varying stiffness k{x) and cross- 
section at rest Aq (x). By a combination of the mass and momentum equations 
in (1), under some regularity assumptions, we have: 

a,. + a. (y + (■'> 

with Ao (x) = k (x) \JAq (x). Considering a stationary flow where the viscous 
friction is neglected by setting Cf = 0, we recover Bernoulli’s law (13). The 
notable difference is that k is now a function of x. In the case of the ’’man at 
rest” equilibrium” (without artifacts such as [26, 35]) we obtain: 

j Q = u = 0 

j dpp (k (x) VA - Ao (x)^ = 0 . 

The fact that now fc is a function of x will influence the way the well-balanced 
scheme is obtained. In the following section, we will present a well-balanced 
scheme for system (1), based on the hydrostatic reconstruction for Saint-Venant/shallow 
water equations with variable pressure [5]. 


4 The numerical method 


4.1 Numerical context 


Several numerical methods have been used to solve the blood flow equations. 
In [43], they are solved thanks to the Methods of Characteristics (MOC). In 
[51, 50], they use a conservative form of the model 

( dtA + dx{Au) = 0 
P 

with the non-conserved vector (A, u) and equations (19) are solved with a two- 
step Lax-Wendroff scheme. In [42], a quasi conservative form of the equations 
(with s{U) a source term) 


= -Cf 


(I 

A2 


(19) 


dtu + dx 


dtA + dxQ — 0 

dtQ + dx = s{U) , 


( 20 ) 
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is solved thanks to a first order explicit in time upwind finite difference scheme. 
In [38], they are the first to solve blood flow equations under a conservative 
form, thanks to a two-step Lax-Wendroff scheme. The solutions of the equations 
under the form (19) using an upwind Discontinuous Galerkin method (used by 
[49, 48]) and a Taylor Galerkin finite element method (also used in [33, 14, 34]) 
have been compared in [41] . A MacCormack finite difference method has been 
applied in [15] followed by [40]. Finite volume methods seem to be first used to 
solve these equations in [9, 10]. In [13], a well-balanced finite volume method 
based on the hydrostatic reconstruction (introduced in [1]) is applied on system 
(5), and this method is compared with a Taylor Galerkin method in [46]. We 
will present in the following sections the extension of the well-balanced scheme 
(based on an extension of the hydrostatic reconstruction) we have used to solve 
the system (1), which can be written under the following vectorial form 


dtU + d,F{U,Z) = Si{U,Z) + S2{U), 


( 21 ) 


with 


U = 


A 


Z = 


Q J ’ 

and the source terms 


Aq 

k 


Q 

, FiU,k)=\ ^ ^^^3/2 

A ^y/np 


( 22 ) 




and S 2 {U) = 


-C 


Q 

A 


(23) 


4.2 Convective step 

For the homogeneous system 


dtU + d,F(u,z) = 0 


(24) 


which is (21) without source term, an explicit first order in time conservative 
scheme can be written as: 


(25) 


where i refers to the cell C = (xi- 1 / 2 , Xi+ 1 / 2 ) = (a;i_i/ 2 , Xi-ij 2 + Ax) and n to 
time tn with tn+i — tn = At. G" is an approximation of U: 




fXi+1/2 


U{x, tn)dx . 


Xi-1/2 


and F^_^i is an approximation of the flux function F(U, Z) at the cell interface 
i + l/2 ' 

F;'Vi/2 = F(c/r,c/r+i,z„z,+i). 

This numerical flux will be detailed in subsection 4.4. 
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4.3 Source terms treatment 

4.3.1 Topography source term Si {U,Z) 

In the system (21), the term Si {U, Z) is involved in the steady state preserva¬ 
tion, therefore requires a well-balanced treatment. Following a variant of the 
hydrostatic reconstruction [5, p.93-94], the variables are reconstructed locally 
from (18) on both sides of the interface i -I- 1/2 of the cell C^: 

{ i/A+i/ 2 L = max(fcV4i -f min(AAoi+i/2,0),0)/fc*^^/2 
Ui+l/2L = (4i+i/2L, 4j_|_i/2L-'Ui)* 

^/A+i/ 2R = max{ki+ii/A^ - max{AAoi+l/2,0),0)/k*_^_J^/^ 

Ui+l/2R = (4i+i/2K, 4i_|_i/2i?-Wi+l)* , 

with AAoi+i /2 = Aoi+i-Aoi = ki+i^/Ao^-k^y/^i and = max(A;i, ki+i). 

In order to help the understanding of the principle of the hydrostatic recon¬ 
struction (26), we present the hydrostatic reconstruction for the shallow water 
system of equations (6): 


hi+i/2L = max(/ii + Zi- z^+ 1 / 2 , 0) 
Ui+i/2L = (^i-|-l/2L) ^i+l/2L-'«i)* 

hi+i/2R = max(/ii+i -I- Zi+i - Zi+1/2, 0) 

Ui+l/2R = ihi+i/2R,hi^i/2R-Ui^iy , 


(27) 


with 2 ;i_|_i /2 = R^3x{zi, Zi+i). The water height is reconstructed in a way that 
allows to have locally the hydrostatic equilibrium h + z = cst on each side of the 
interface j -I- 1/2. As mentioned in [1], max(., 0) is there to ensure the positivity 
of the water height in case of drying and the upwind evaluation of 2 i+i /2 ensures 
that 0 < hi+i/ 2 L A hi and 0 < /ii+i/ 2 i? < ^j+ij which has been proved in [1] 
to ensures the positivity of the water height. For blood flow equations with a 
constant stiffness k, the corresponding equilibrium writes -s/A — = cst, so 

\fA (respectively —^/Ao) ’’plays the role” of h (resp. z), thus in that case the 
hydrostatic reconstruction writes: 



(28) 


As we have instead of z, we take = min(-v/Ao~, i/A qi^i), thus 

we have: 


\/Ai+i/2L = max(VA/ -h min(AVA(j+i/2,0), 0) 

Ui+ll2L = (Ai-|-i/2L) Ai-|-i/ 2Z,-'^i)* 

i/Ai+i/2R = max(y^Ai+i - max(A^/Aoj+i/2, 0), 0) 

Ui+l/2R = (Ai+1/2R, Aj_|_i/2fi-'Wi+l)‘ , 


(29) 


with Ay/Aiiij^ii 2 = \/Aoj+i — y/Aoi. We can notice that we recover reconstruc¬ 
tion (29) if the stiffness k is constant in reconstruction (26). For consistency, 
the scheme (25) is modified as follows: 



( 30 ) 


















where 



with 




Si+l/2L 



and P(A, fc) = k{x) /{ipy/Tr). Thus blood flow in a artery with varying 
cross-section at rest and stiffness is treated in a well-balanced way. 

4.3.2 Viscous source term S 2 (U) 

In system (21), the friction term —CfQjA in S '2 (U) is treated semi-implicitly. 
This treatment is classical in shallow water simulations [7, 31] and has proven 
efficient in blood flow simulation as well [13]. Furthermore, this treatment 
preserves the ’’dead man” equilibrium (18). It consists in using first (30) as a 
prediction step without friction, i.e.: 



then applying a semi-implicit friction correction on the predicted values (U*): 



Thus we get the corrected velocity and we have A^^^ = A*. 

4.4 HLL numerical flux 

As presented in [13], several numerical fluxes can be used (Rusanov, HLL, 
VFRoe-ncv and kinetic fluxes) for numerical simulations of blood flow in ar¬ 
teries. Details can be found in [5, 12, 13]. In this work we will use the HLL 
flux (Harten Lax and van Leer [21]) because it is the best compromise between 
accuracy and CPU time consumption (see [12, chapter 2]). It writes: 


F(UL,t/fl,fc*) 


F{UL,k*) 

C2F{UL,k*)-ciFiUR,k*) ^ 


C1C2 


{Ur - Ul) if Cl < 0 < C2 
if C2 < 0 , 


if 0 < Cl 


C2 - Cl 


C2 - Cl 


F{UR,k*) 


with 


Cl = inf ( inf A, (C/, A:*)) and C 2 = sup ( sup Xj{U,k*)), 



where Xi{U,k*) and X 2 {U,k*) are the eigenvalues of the system and k* 
max(fcL, kn). 
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To prevent a blow up of the numerical values, we impose the following CFL 
(Courant, Friedrichs, Lewy) condition: 


At < ncFL 


Ax 


max(|ui| + Ci) ’ 


where c* = ^kiy^i/{2py^) and ucfl = 1- 


5 Validation of the method 

To validate the well-balanced scheme presented in the previous sections for blood 
flow in arteries with varying stiffness k (x) and cross-section at rest (3^)? we 
applied it to different test cases taken from [13], where arteries with a varying 
cross-section at rest Aq (x) and a constant stiffness k were considered. For each 
of these examples, the rest equilibrium state was: <5 = 0 and \fA — y/A^ = 0 
and non-reflecting boundary conditions were set at each end of the computa¬ 
tional domain in the form of homogeneous Neumann boundary conditions. The 
hydrostatic reconstruction scheme as well as a naive centered discretization of 
the source term were systematically tested to clearly evaluate the benefit of us¬ 
ing a well-balanced scheme. According to [13], several Riemann solvers can be 
used, but we only display results obtained using the HLL flux. In the following, 
we present the numerical parameters, the analytic solution if it exists and the 
numerical results. For further details we refer the reader to [13]. 


5.1 ’’The man at eternal rest” 

We considered an artery at its equilibrium state, where there is no flow and 
the radius of the cross-section at rest Rq {x) varies throughout the artery, as for 
example in a dead man with an aneurysm. This equilibrium state is exactly the 
one well-balanced methods are designed to preserve. If the topography source 
term is not treated correctly, non-physical velocity may be generated. 

We used the following numerical values: L = 0.14m,J = 50 cells. Tend = 
5 s, p = 1060 kg.m~^, C/ = 0 and fc = 4.0 x 10® Pa.m~^. We used the 
equilibrium state as an initial condition, setting Q(x,0) = 0 and: 


R {x, 0) = Rq (x) = < 


Rq 
Rq 
Rq 
Rq 

. Rq 

1-3 


AR 

AR 

AR 


1 -I- sin ( - - -I- TT 


X — Xi 
X2 - Xi 


1 -I- cos TT 


X - Xq 
X i - Xq 


if a; S [0, Xi] 
if a; e ]a;i, a: 2 [ 
if a: G [a: 2 ,a: 3 ] 
if a: G ]a:3, x^l 
if a; G [xi, L] , 


with Rq = 4.0 X 10“® m, AR = 1.0 x 10“® m, a:i = 1.0 x 10“^ m, x^ = 
3.05 X 10“^ m, 3:3 = 4.95 x 10“^ m and x^ = 7.0 x 10“^ m. The radius at rest 
is plotted on Figure 2 left. 
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The results obtained are presented in Figure 2 right. As expected, a naive 
centered discretization of the topography source term results in nonphysical os¬ 
cillations of the velocity u(x,t), whereas the well-balanced solution preserves 
the equilibrium state. 



1.5 
1 

0.5 
0 

■0.5 
-1 
■1.5 

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 

x[m] 


A Hydrostatic reconstruction - 

/* Centered discretization ■■■■ 

t • 

« • 
a a 
a a 



Figure 2: The ’’dead man case”: (Left) The radius of the artery Rq {x); (Right) 
Comparison of the velocity at time t = 5 s between an explicit treatment of the 
source term (dashed line) and the hydrostatic reconstruction (full line). 


5.2 The ideal ’’Tourniquet” 

This test case is the equivalent of the dam break problem for the Shallow Water 
equations (Stoker’s solution in [12]). We considered an artery with a constant 
radius at rest Rq, a constant stiffness k and no viscous friction (C/ = 0), there¬ 
fore the governing system of equations was (24). Initially, a tourniquet was 
applied and then immediately removed. We have a Riemann problem and the 
method of characteristics allowed us to compute an analytic solution that we 
compared to the numerical solutions. This Riemann problem has been first in¬ 
troduced in compressible gas dynamic with the Sod tube (for further details we 
refer the reader to [28, 32]) and extended to blood flow in [13]. 


We considered an artery of length L = 8.0 x 10“^ m with x G > f] and 
used the following numerical parameters: J = 100 cells. Tend = 5.0 x 10“^ s, 
p = 1060 kg.m~^ and k = 1.0 x 10^ Pa.m~^. We used a perturbation of the 
equilibrium state as an initial condition, setting Q(x,0) = 0 and: 


A (x, 0) 


Al =7r {Rq + AR)^ 
An =7ri?o 


if a; G 
if a; G 


“■f 


with Rq = 4.0 X 10 ^ TO and AR = 1.0 x 10 ^ to. 
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The results obtained are presented in Figure 3. We can see that the nu¬ 
merical solution obtained with the well balanced scheme is in good agreement 
with the analytic solution presented in [13]. This is also true for the solution 
obtained using a centered discretization of the topography source term, which 
is superposed on the well-balanced solution, since in this case the source term 
is null. 



X [m] X [m] 


Figure 3: The Tourniquet: (Left) Radius of the artery R{x) at t = 5 x 10“^ s; 
(Right) Flow rate of the artery Q{x) at t = 5 x 10“^s. Comparison between 
the exact analytic solution (full line) and the numerical solution obtained with 
an explicit treatment of the topography source term and the hydrostatic recon¬ 
struction (dashed lines). The numerical solutions are superposed. 


5.3 Wave reflection-transmission of the pulse towards a 
constriction 

In this section we considered the propagation of a pulse towards constriction. 
This configuration is an idealized representation of a transition between a parent 
artery and a daughter artery of smaller cross-section. We tested here the ability 
of the numerical scheme to capture the propagation of a small perturbation of 
the equilibrium state at the beginning of an artery with a varying radius at 
rest Ro{x). In order to accurately compute the numerical solution, the forward 
and backward traveling waves need to be correctly captured as well as the re¬ 
flected and transmitted waves generated by the abrupt change in topography at 
the transition point. To test if these reflections were accurately described, we 
computed the analytic reflection and transmission coefficients at the transition 
point and compared them to the amplitude of the numerical reflected waves. 
For further details we refer the reader to [13]. 

We considered an artery of length L = 0.16 to and used the following numer¬ 
ical parameters: J = 1500 cells. Tend = 8.0 x 10“^ s, p = 1060 kg.m~^, C/ = 0 
and fc = 1.0 X 10® Pa.m~^. The constriction was defined by the following radius 
of the cross-section at rest: 
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Rq (x) 


Rn + Ai? 
= {Rr + ^ 


1 + COS TT 


X — XI 
X2 - Xi 


if X S [ 0 , Xi] 
if X e ]xi,X 2 ] 
if X G ]x2, L] , 


with Rr = 4.0 X 10 ^ TO, Ai? = 1.0 x 10 ^ to, xi = and X 2 = f. We set 
Q(x, 0 ) = 0 as an initial condition and we defined the initial perturbation as: 


R (x, 0) 


I Ro{x) 


^100 , A' 

1 - 1 - e sin 


[ Ro{x) 




if X G [x 3 , X 4 ] 

else , 


with X 3 = < xi, X 4 = < X 2 and e = 5.0 x 10“^ a small parameter 

ensuring that we stayed in the range of small perturbations of the equilibrium 
state. 


The numerical results are plotted in Figure 4. We can see that the propaga¬ 
tion of the pulse as well as the wave reflections and transmissions are accurately 
described using the well balanced scheme (Figure 4 left) whereas spurious waves 
appear with the centered discretization of the source term (Figure 4 right). 



X [m] X [m] 

Figure 4: (Left) Hydrostatic reconstruction; (Right) Centered discretization of 
the topography source term. R{x) — Rq{x) at 3 time steps: t = 0, t = t = 
The straight dashed lines represent the level of the predicted reflection 
{Re) and transmission (T^) coefficients. 


6 Asymptotic solutions for a uniform vessel 

In this section we studied the propagation of a pulse wave in a uniform vessel 
{k = cst, Aq = cst) and derived asymptotic solutions of the system of equation 
(1), following the work of Wang and al. [47]. Small perturbations ( eQ, Aq + eAj 
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of the base state {Q = Q, A = Ag) were considered, resulting in the following 
linearized system of equations: 



(31) 


where cg = ^/kRo/ (2p) is the Moens-Korteweg celerity. 

In the following numerical examples, we only present results obtained for 
the hydrostatic reconstruction since we considered a uniform vessel. The nu¬ 
merical parameters were defined as follows: L = 3 m, Rq = 1.0 x 10“^ m, 
J = 1500 cells. Tend = 0.5 s, p = 1060 kg.m~^, p = 3.5 x 10“^ Pa.s and 
/c = 1.0 X lO’^ Pa.m~^. The parameters C/ and Cy, respectively the viscous 
coefficient and the viscoelastic coefficient, were set according to the desired test 
case. 

Initially, the system was at its equilibrium state (Q = 0, A = Ag = ttRq) and 
an inflow boundary condition was prescribed as Q = 0, t) = Qin it) with: 



where Hit) is the Heaviside function, Te the period of the sinusoidal wave and 
Qc the maximum amplitude of the inflow wave. We set Qc = 1-0 x 10“® m?.s~^ 
and Tc = 0.4 s to insure that only small perturbations from the equilibrium state 
were considered. The cross-section at the inlet A{x = 0, t) was reconstructed by 
a matching of the outgoing characteristic, technique that takes advantage of the 
hyperbolic nature of the problem. A homogeneous Neumann boundary condi¬ 
tion was prescribed at the outlet to simplify the computation of the asymptotic 
solutions and to avoid reflections. 

6.1 The d’Alembert equation 

Following ideas developed in [47], we set (7/ = 0 in (31) and we obtained the 
d’Alembert equation, which admits the following pure wave solution cgAg = 
Q — Qin Cgt). 

In Figure 5, we can see the propagation of a pulse wave without dissipation 
or diffusion, as predicted by the analytic solution. 

6.2 Dissipation due to the viscosity of the blood 

We also investigated the effect of the blood viscosity on the propagation of the 
pulse wave and set Cf ^ 0. Starting from the linearized system of equations 
(31), we considered the small parameter £/ = and performed the change 
of variables ^ = x — c^t and t = eft to place ourselves in the moving frame at 
slow times to properly capture the effects of the viscous term. The first order 
solution obtained in [47] is: 
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Figure 5: Pure wave solution: u(x) at time t = {0.1, 0.2,0.3, 0.4,0.5} for the 
well-balanced scheme. The straight black dotted line represents the maximum 
amplitude of the pure wave solution. 


cqAo = Qo = Qo (x - Cot) exp 

where exp is the exponential envelop of the pure wave solution Qq (x — cgt). 

To obtain this asymptotic solution numerically, we set Cf = AOtti/ = 4.15 x 
10“^ therefore e/ = 0.53. 

In Figure 6, we can see the propagation of the pulse with dissipation (or 
attenuation) of its amplitude due to the viscosity of the blood. The straight 

doted line represents the exponential envelop exp computed previ¬ 

ously and is in good agreement with the decrease in amplitude of the pulse wave. 

One can note that as expected, there is no diffusion, since the wavelength of the 
pulse does not change while it propagates in the artery. 


6.3 Diffusion due to the viscoelasticity of the arterial wall 

In this section, we set the friction coefficient to zero {Cf =0) and focused on an 
other important characteristic of the blood flow in the arteries: the viscoelastic¬ 
ity of the arterial wall. We chose here to take into account this time-dependent 
behavior in our governing system of equations through a very simple lumped 
model, the Kelvin-Voigt model, resulting in an additional parabolic term in the 
governing system of equations: 


dtA + dxQ — 0 

dtQ + d,(^ + ’ 


(32) 


where the viscoelastic coefficient is defined as = 1-57 rri^-s ^ 

with (j) = 5000 Pa.s and h = 5.0 x 10“^ m. The parabolic term was treated 
by performing a temporal splitting of the problem. First the purely hyperbolic 
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Figure 6: Viscous damping: u{x) at time t = {0.1,0.2, 0.3,0.4, 0.5} for the 
well-balanced scheme. The straight black dotted line represents the exponential 
envelop of the asymptotic solution. 


problem with a non reflecting boundary condition at the outlet was solved, and 
its solution was then used as an initial condition of the parabolic problem. A 
Crank-Nicolson scheme coupled with homogeneous Neumann boundary condi¬ 
tions was than used to solve the parabolic problem. 


To correctly capture the behavior of this new viscoelastic term, we defined 
a new small parameter Ci, = = 8.3 x 10“^ and applied the same technique 

CqIc 

as in the previous section. From [47] we have the following first order diffusive 
analytic solution, which is a solution of the heat equation: 


[ Qo(t,C) 

[G(r,e) = 



Qo (0, Jy) G (t, ^ - ry) dry 


^ p-CV(2TegT„) ^ 

^J2tttcITc 


The numerical results for several times and the analytic solution at t = 0.4 s 
are presented in Figure 7. We can see that the viscoelastic term induces a 
diffusion of the pulse wave, changing its wavelength, and that the numerical 
solution at t = 0.4 s perfectly matches with the asymptotic solution at t = 0.4 s. 


7 Real artery simulation 

In this section, we focused on simulating the propagation of a pulse wave in a 
tapered artery of length L = 3 m, where the the radius of the cross section at 
rest Ro{x) was linearly decreasing from the proximal to the distal end of the 
artery: 


r rl 

Ro{x) = < Rl — [x — xi)AR 
[ - {X2 - Xi)AR 
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if x G [0, a;i[ 
if x G [xi,X2[ 
if X G [x2, L[ , 
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Figure 7: Viscoelastic diffusion; u{x) at time t = {0.1,0.2,0.3,0.4,0.5} for 
the well-balanced scheme (dashed lines). The black dotted line represents the 
asymptotic solution at t = 0.4 s. 


with Rl = 4.0 X 10“^ m, AR = 1.0 x 10“^ m, Xi = and X 2 = Follow¬ 
ing [47], the stiffness of the arterial wall was defined as k{x) = | with E the 
Young’s modulus and h the width of the arterial wall. Therefore we were in a 
configuration where i?o and k were varying throughout the length of the artery 
and if the well-balanced scheme was not used, spurious waves might have arisen. 

We used the following numerical parameters to mimic the geometrical and 
mechanical properties of a real artery: J = 1500 cells. Tend = 0.5 s, p = 
1060 kg.m~^, p = 3.5 x 10“^ Pa.s, E = 4.0 x 10® Pa, h = 5.0 x lO”’^ to, 
Cf = Sttv, (p = 5000 Pa.s and We used the same initial inflow 

condition as for the asymptotic solutions. 



X [m] X [m] 


Figure 8: Tapered artery - Pure wave solution: u {x) at time t = 
10.1,0.2,0.3,0.4,0.5} for Cf = 0 and = 0: (Left) Centered discretization 
of the topography source term; (Right) Hydrostatic reconstruction. 
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Figure 9: Tapered artery: viscous and viscoelastic effects: u {x) at time t = 
{0.1, 0.2,0.3, 0.4,0.5} for the well-balanced scheme. 


The results are presented in figures 8 and 9. We can see that in the absence 
of friction and viscoelastic effects (figure 8), if the well-balanced scheme is not 
used (figure 8 left) nonphysical reflections appear. On the contrary, the well- 
balanced scheme provides a satisfactory numerical solution, where a continuous 
reflection phenomena takes place due to the tapering, resulting in a decrease of 
the amplitude of the backward traveling wave and an increase of the amplitude 
of the forward traveling wave. Indeed, in the case of a tapered artery, the 
transmission coefficient > 1 and the reflection coefficient i?e < 1. When 
viscous and viscoelastic effects are taken into account (figure 9), all phenomena 
add up and we recognize the effects of the continuous reflection, the viscous 
dissipation and the viscoelastic diffusion. 


Conclusion and perspectives 

In this work we have presented a numerical method based on a well-balanced 
finite volume scheme for the blood flow equations with variable wall elasticity. 
This scheme based on an extension of the hydrostatic reconstruction gave very 
good results on several tests, for which classical methods failed. In further work, 
we will try to improve the accuracy of the numerical method by raising the order 
of the numerical method and to apply this method to real network modeling. 
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